function [label, markslbl, marks, obj, Zv, alpha] = MVSC(data, k, opts)
% [label, obj] = MVSC(data,k,opts)
% 
% Mutli-view Spectral Clustering
% 
% Input:
%       - data: nv cell aray of the data matrix of size nSmp x nFea, 
%               where each row is a sample point
%       - k: the number of clusters
%       opts: options for this algorithm
%           - p: the number of salient points picked (default 1000)
%           - r: the number of nearest salient points for representation (default 5)
%           - wr: the exponential of the weight between views (default 1)
%           - numRep: the number of replicates for the final kmeans (default 10)
%           - maxIter: the maximum number of iterations for final kmeans (default 100)
%           - mode: landmark selection method, currently support
%               - 'kmeans': use centers of clusters generated by kmeans (default)
%               - 'random': use randomly sampled points from the original
%                           data set 
%           The following parameters are effective ONLY in mode 'kmeans'
%           - kmNumRep: the number of replicates for initial kmeans (default 1)
%           - kmMaxIter: the maximum number of iterations for initial kmeans (default 5)
% Output:
%       - label: the cluster assignment for each point
% Requre:
%       litekmeans.m
% Usage:
%       data = {rand([100,50])};
%       label = MVSC(data,10);
%
%   Written by Yeqing Li (yeqing.li AT mavs.uta.edu)
%              
% Set and parse parameters
if (~exist('opts','var'))
   opts = [];
end


p = 1000;
if isfield(opts,'p')
    p = opts.p;
end

r = 5;
if isfield(opts,'r')
    r = opts.r;
end

maxIter = 100;
if isfield(opts,'maxIter')
    maxIter = opts.maxIter;
end

numRep = 10;
if isfield(opts,'numRep')
    numRep = opts.numRep;
end

mode = 'kmeans';
if isfield(opts,'mode')
    mode = opts.mode;
end

wr = 1;
if isfield(opts,'wr'),
    wr = opts.wr;
end

maxWghtIter = 5;
if isfield(opts,'maxWghtIter'),
    maxWghtIter = opts.maxWghtIter;
end

thresh = 1e-6;
if isfield(opts,'thresh'),
    thresh = opts.thresh;
end

kertype = 'Gaussian';
if isfield(opts,'kertype'),
    kertype = opts.kertype;
end

%===================================================================

nSmp=size(data{1},1);
numView = length(data);
% inti alpha
alpha = ones(1, 1, numView)/numView; 

%===================================================================
% Salient point selection
if strcmp(mode,'kmeans')
    conD = cell2mat(data);
    nCols = cell2mat(cellfun(@(x) size(x, 2), data, 'UniformOutput', 0));
    kmMaxIter = 5;
    if isfield(opts,'kmMaxIter')
        kmMaxIter = opts.kmMaxIter;
    end
    kmNumRep = 1;
    if isfield(opts,'kmNumRep')
        kmNumRep = opts.kmNumRep;
    end
    [dump,marks]=litekmeans(conD,p,'MaxIter',kmMaxIter,'Replicates',kmNumRep);
    marks = mat2cell(marks, p, nCols);
    clear kmMaxIter kmNumRep
elseif strcmp(mode,'random')
    indSmp = randperm(nSmp);
    marks = cell(1, numView);
    for i = 1:numView,
        marks{i} = data{i}(indSmp(1:p),:);
    end
    clear indSmp
else
    error('mode does not support!');
end

%===================================================================
% Z construction for multi-view
Zv = zeros(nSmp, p, numView);
for iView = 1:numView,
    if strcmp(kertype, 'Gaussian'),
        D = EuDist2(data{iView},marks{iView},0);

        if isfield(opts,'sigma')
            sigma = opts.sigma;
        else
            sigma = mean(mean(D));
        end

        dump = zeros(nSmp,r);
        idx = dump;
        for i = 1:r
            [dump(:,i),idx(:,i)] = min(D,[],2);
            temp = (idx(:,i)-1)*nSmp+[1:nSmp]';
            D(temp) = 1e100; 
        end

        dump = exp(-dump/(2*sigma^2));
    elseif strcmp(kertype, 'Linear'),
        sigma = optSigma(data{iView});
        options.KernelType = 'Linear';
        options.t = sigma; % width parameter for Gaussian kernel
        Sim = constructKernel(data{iView},marks{iView},options);
        
        dump = zeros(nSmp,r);
        idx = dump;
        for i = 1:r
            [dump(:,i),idx(:,i)] = max(Sim,[],2);
            tep = (idx(:,i)-1)*nSmp+[1:nSmp]';
            Sim(tep) = 1e-100;
        end
        %ind = (idx-1)*n+repmat([1:n]',1,s);
    else
        error('Unknown kernel type');
    end
    
    sumD = sum(dump,2);
    Gsdx = bsxfun(@rdivide,dump,sumD);
    Gidx = repmat([1:nSmp]',1,r);
    Gjdx = idx;
    Z = sparse(Gidx(:),Gjdx(:),Gsdx(:),nSmp,p);
    
    feaSum = full(sqrt(sum(Z,1)));
    feaSum = max(feaSum, 1e-12);
    Z = Z./feaSum(ones(size(Z,1),1),:);
    
    Z(isinf(Z)) = 0;
    Z(isnan(Z)) = 0;
    
    Zv(:, :, iView) = Z;
end
% Z = sum(Zv, 3);
tmp = 1/(1-wr);
obj = zeros(maxWghtIter, 1);
for t = 1:maxWghtIter,
    t
    Z = sum(bsxfun(@times, Zv, alpha.^wr), 3);
    Z = sparse(Z);

    %===================================================================
    % Graph decomposition
    [U, S, V] = mySVD(Z,k+1);
    U(:,1) = [];
    V(:,1) = [];

    U=U./repmat(sqrt(sum(U.^2,2)),1,k);
    V=V./repmat(sqrt(sum(V.^2,2)),1,k);
    
    % Update alpha
    if ~isinf(tmp) && ~isnan(tmp)
        for j = 1:numView,
            h(j) = trace(U'*Zv(:, :, j)*Zv(:, :, j)'*U);
            %h(j) = trace(U'*Zv(:, :, j)*Zv(:, :, j)'*U) + trace(V'*Zv(:, :, j)'*Zv(:, :, j)*V);
        end
    
        alpha = ((wr*h).^tmp)/(sum(((wr*h).^tmp)));
        alpha = reshape(alpha, [1 1 numView]);
    end
    
    % calculate the obj
    obj(t) = 0;
    for v = 1: numView
        obj(t) = obj(t) + (alpha(v)^wr)*h(v);
    end
    if(t > 1)
        diff = obj(t-1) - obj(t);
        if(diff < thresh)
            break;
        end
    end
end

% Final kmeans
alllabel = litekmeans([U; V],k,'MaxIter',maxIter,'Replicates',numRep);
label = alllabel(1:nSmp);
markslbl = alllabel(nSmp+1:end);
